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We develop a linear method for solving the nonlinear differential equations of a lumped-parameter 
thermal model of a spacecraft moving in a closed orbit. Our method, based on perturbation the- 
ory, is compared with heuristic linearizations of the same equations. The essential feature of the 
linear approach is that it provides a decomposition in thermal modes, like the decomposition of 
mechanical vibrations in normal modes. The stationary periodic solution of the linear equations can 
be alternately expressed as an explicit integral or as a Fourier series. We apply our method to a 
minimal thermal model of a satellite with ten isothermal parts (nodes) and we compare the method 
04 . with direct numerical integration of the nonlinear equations. We briefly study the computational 

complexity of our method for general thermal models of orbiting spacecraft and conclude that it is 
■ certainly useful for reduced models and conceptual design but it can also be more efficient than the 

<n: direct integration of the equations for large models. The results of the Fourier series computations 

^ . for the ten-node satellite model show that the periodic solution at the second perturbative order is 

^ ■ sufficiently accurate. 
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Ai = outward-facing area of ith-node, re? 

C = thermal capacitance diagonal matrix, J/K 

• Ci = thermal capacitance of zth-node, J/K 
= eigenvector of C^^^JC'^^^, K 
= jth-component of eigenvector Ca, K 

F{t) — driving function vector, K/s 

■ F{m) ~ mth Fomder coefficient of F{t), K/s 
. F* — complex conjugate of F, K/s 

■ Fa{t) = driving function for ath-mode, K/s 
', Fa{t) = time derivative of Fa{t), K/s^ 

■ Flit) = driving function for jth-node, K/s 

. Fi{m) = mth Fourier coefficient of Fi{t), K/s 

" G{t) — driving function for second order temperature vector, K/s 

' , I [ J — Jacobian matrix, s~^ 

^ ■ Jij — Jacobian matrix ij-element, s~^ 

k ~ number of steps in numerical integration of ODEs 

K — conduction coupling matrix, W/K 

" — matrix of conductances from linearized radiation coupling terms, W/K 

— ith-node conductance from linearized environment-radiation terms, W/K 

Kij — conduction coupling matrix ij -element, W/K 

^ij — jj-element of conductance matrix from linearized radiation terms, W/K 

n = number of samples for discrete Fourier transform 

N = number of nodes 

P ~ eigenvector matrix for J, K 

Pia — ith-component of J-eigenvector for ath-mode, K 

q = auxiliary heat input vector, K/s 

Qi{t) = heat input to ith-node, W 

{Qi) — mean value of Qi{t) over period T, W 

Ri = zth-node coefficient of radiation to the environment, W/K* 

Rij — radiation coupling matrix, W/K"* 
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cosmic microwave background radiation temperature, K 
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nth order term in jth-node temperature expansion, K 
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stationary solution temperature vector, K 
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nth order of stationary solution temperature vector, K 
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time derivative of ith-nodc temperature, K/s 
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iiifx'arcd cixiissivity of ztlx-node outward-facing surface 
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perturbation parameter 


Aa 




eigenvalue of J for ath mode, s~^ 


SXa 




perturbation of J-cigcnvaluc for ath mode, s"""^ 


a 




Stefan-Boltzmann constant, 5.67 x 10~^ Wm~^K~^ 



I. INTRODUCTION 

The thermal control of a spacecraft ensures that the temperatures of its various parts are kept within their appro- 
priate ranges [1-4] . The simulation and prediction of temperatures in a spacecraft during a mission are usually carried 
out by commercial software packages. These software packages employ "lumped parameter" models that describe the 
spacecraft as a discrete network of nodes, with one energy-balance equation per node. The equations for the thermal 
state evolution are coupled nonlinear first-order differential equations, which can be integrated numerically. Given the 
thermal parameters of the model and its initial thermal state, the nunicirical integration of the differential equations 
yields the solution of the problem, namely, the evolution of the node temperatures. However, a detailed model with 
many nodes is difficult to handle, and its integration for a sufficiently long time of evolution can take considerable 
computer time and resources. Therefore, it is very useful to study simplified models and approximate methods of 
integrating the differential equations. 

Many spacecraft missions, in particular, satellite missions, consist of an initial transient part and then a stationary 
part, in which the spacecraft just goes around a closed orbit, in which the heat inputs are periodic. These periodic 
heat inputs arc expected to induce periodic temperature variations, with a maximum and a minimum temperature 
in (;ach orbit. This suggests a conservative approach that consists in computing only the temperatures for the hot 
and cold cases of the given orbit, defining them as the two steady cases with the maximum and minimum heat loads, 
respectively. Naturally, the real temperature variations in the orbit are smaller, because there is not enough time for 
the hot and cold cases to establish themselves. In fact, the temperature variations can be considerably smaller, to 
such a degree that it is necessary to integrate the differential equations, at least approximately. 

The diff'erential equations for energy balance are nonlinear due to the presence of radiation couplings, which follow 
the Stefan-Boltzmann quartic law. A common approach to these equations involves a linearization of the radiation 
terms that approximate them by heat conduction terms [3, 5-7]. This approach transforms the nonlinear equations 
into standard linear heat conduction equations. But this approach has not been sufficiently justified, is of a heuristic 
nature and does not constitute a systematic approximation. 

In fact, nonlinear equations are very different from linear equations and, in particular, a periodic driving may not 
induce periodic solutions but much more complex solutions, namely, chaotic solutions. Therefore, we have carried out 
in preceding papers a full nonlinear analysis of spacecraft thermal models [8, 9]. The conclusion of the analysis is that 
the complexities of nonlinear dynamics, such as multiple equilibria and chaos, do not appear in these models. While 
the existence of only one cqiiilibrium state can be proved in general, the absence of chaos under driving by variable 
external heat loads can only be proved for a limited range of magnitudes of the driving loads. This range presumably 
includes the magnitudes involved in typical spacecraft orbits. The proofs in Refs. 8 and 9 are constructive and are 
based on a perturbation method that is expected to be sound when the linear equations corresponding to the first 
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perturbativc order constitute a good approximation of the nonlinear equations. This imphes that the fuhy nonhnear 
solution describes a weakly nonlinear oscillator. Since the perturbative approximation is mathematically rigorous and 
systematic, it is worthwhile to study in detail the scope of the perturbative linear equations and, furthermore, to 
compare them with previous linear approaches of a heuristic nature. 

The main purpose of this paper is to study the linear method of predicting the thermal behavior of spacecraft in 
stationary orbits (Sect. II and III) and to test it on a minimally realistic thermal model of a satellite in a circular 
orbit. Since the general one and two-node models analyzed in Refs. 8 and 9, respectively, are too simple, we define 
in this paper a ten-node thermal model of a small Moon-orbiting satellite (Sect. IV). This model is simple enough 
to allow us to explicitly show all the quantities involved (thermal couplings and capacities, heat inputs, etc.) and it 
is sufRcient for illustrating the main features of the linear approach. As realistic thermal models have many more 
nodes, we consider in Sect. V the important issue of scalability of the method and, hence, its practical applications. 
Computational aspects of the steady-state problem have been studied by Krishnaprakas [10, 11] and by Milman and 
Petrick [12], while computational aspects of the direct integration of the nonlinear equations for the unsteady problem 
have been studied by Krishnaprakas [13]. Here we focus on the linear equations for the stationary but unsteady case 
and survey its computational aspects. 

A note on notation: In the equations that contain matrix or vector quantities, sometimes we use component notation 
(with indices) while other times we use compact matrix notation (without indices), according to the nature of the 
equations. 

II. LINEARIZATION OF THE HEAT-BALANCE EQUATIONS 

A lumped-parameter thermal model of a continuous system consists of a discrete network of isothermal regions 
{nodes) that represent a partition of the total thermal capacitance and that are linked by thermal conduction and 
radiation couplings [1-5]. This discretization reduces the integro-differential heat-transfer equations to a set of energy- 
balance ODEs, one per node, which control the evolution of the nodes' temperatures [5]: 

N 

Ci ti = Qi{t) - ^ [Kij{Ti - Tj) + Rij{T^^ - T/)] - Ri {Ti^ - T^), i = l,...,N, (1) 
i=i 

where N is the number of nodes and Qi{t) contains the total heat input to the zth-node from external radiation and 
from internal energy dissipation (if there is any). The conduction and radiation coupling matrices are denoted by 
K and i?, respectively; they are symmetric {Kij — Kji and Rij = Rji) and Ku = Ru = 0; so there are A'^(A'^ — 1) 
independent coupling coefficients altogether, but many vanish, usually. The temperature T}) ~ 3K is the temperature 
of the environment, namely, the cosmic microwave background radiation. The ith-node coefficient of radiation to the 
environment is given by Ri = AiSia, where Ai denotes the outward facing area, Si its (infrared) emissivity, and a is the 
Stefan-Boltzmann constant. The constant term RiT^ can be included in Qi{t) or ignored altogether, if each Tj ^ Tq. 
Equations (1) coincide with the ones implemented in commercial software packages, for example, ESATAN"""^ [14]. 

There is no systematic procedure for finding the analytical solution of a system of nonlinear differential equations, 
except in some particularly simple cases. Of course, nonlinear systems can always be integrated numerically with 
finite difference schemes. Methods of this kind are employed in commercial software packages. When a nonlinear 
system can be approximated by a linear system and, hence, an approximate analytic solution can be found, this 
solution constitutes a valuable tool. Actually, one can always resort to some kind of perturbation method to linearize 
a nonlinear system. Therefore, we now study the rigorous linearization of Eqs. (1) based on a suitable perturbation 
method, and we also describe, for the sake of a comparison, a heuristic linearization, which actually is best understood 
in light of the results of the perturbation method. 

A. Perturbative lineEirization 

If we assume that the heat inputs Qi{t) in the energy-balance Eqs. (1) are periodic, namely, that there is a time 
interval T such that Qi{t + T) = Qi{t), then it seems sensible to study first the effect of the mean heat inputs in a 
period. This averaging method, introduced in Refs. 8 and 9, relies on the fact that the autonomous nonlinear system of 
ODEs for constant Qi can be thoroughly analyzed with analytical and numerical methods. For example, it is possible 
to determine that there is a unique steady thermal state and that it is (locally) stable [9, 12]. The actual values of the 
steady temperatures can be found efficiently with various numerical methods [10-12]. Furthermore, the eigenvalues 
and eigenvectors of the Jacobian matrix of the nonlinear system of ODEs provides us with useful information about 
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the dynamics, in particular, about the approach to steady-state: the eigenvectors represent independent thermal 
modes and the eigenvalues represent their relaxation times [9]. 

Once the averaged equations arc solved, the variation of the heat inputs can be considered as a driving of the 
averaged solutions. Thus, we can define the driving function 



m) = 



Q^it) - (Qi) 

Ci 



i = l,...,N, 



where (Qi) denotes the mean value of Qi(t) over the period of oscillation. A weak driving function must not produce 
a notable deviation from the averaged dynamics. In particular, the long-term thermal state of an orbiting spacecraft 
must oscillate about the corresponding steady-state. To embody this idea, we introduce a formal perturbation 
parameter e, to be set to the value of unity at the end, and write Eqs. (1) as 



m 

Ci 



N r 



E 



1,...,7V, 



Then, we assume an expansion of the form 



(2) 



Ti{t) = Y,e^T(^)^{t). 



(3) 



When we substitute this expansion into Eqs. (2), we obtain for the zeroth order of e 

N r - 
J 

'Ci 



T, 



(0)z 



Ci 



E 



^i-'(0)i - -t(0)jj + 7^(,J(0)i 



^(0)i) 



Ri 



l,...,Ar, 



(4) 



that is to say, the averaged equations. The initial conditions for these equations are the same as for the unaveraged 
equations. 

For the first order in e, we obtain the following system of linear equations: 



N 



t(l)i=Y,Ji3{t)T(l)0+Fi{t): 



i = l,...,N. 



(5) 



Here, Jij{t) is the Jacobian matrix 



Jij{t) 



Ti{T) 



where T(^o){t) is the solution of the zeroth order equation. Equations (5) are to be solved with the initial condition 

T(i)(0)=0. 



The elements of the Jacobian matrix at a generic point in the temperature space are calculated to be: 

Jij = Cr^ {Kij + 4RijTf) , if i ^ j, 



J a — Cj 



N 



■Y,{Kik+^RikTf)-ARiT;^ 



L fc=l 



(6) 
(7) 



This matrix has interesting properties. First of all, it has negative diagonal and nonnegative off-diagonal elements. In 
other words, — J is a Z-matrix [15]. Furthermore, it fulfills a semipositivity condition that qualifies it as a nonsingular 
Af -matrix [9]. Since the eigenvalues of an M-matrix have positive real parts, the opposite holds for J, namely, 
its eigenvalues have negative real parts. One more interesting property of — J, related to semipositivity, is that it 
possesses a form of diagonal dominance: it is similar to a diagonally dominant matrix and the similarity is given by a 
positive diagonal matrix. Naturally, this property is shared by J. These properties are useful to prove some desirable 
properties of the solutions of Eqs. (5). 

The chief property of J is that — J is a nonsingular M-matrix. In particular, it implies that —J~^ is non- negative 
and, therefore, that the Perron- Frobenius theory is applicable to it [15]. The relevant results to be applied are: (i) 
Perron's theorem, which states that a strictly positive matrix has a unique real and positive eigenvalue with a positive 
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eigenvector and that this eigenvalue has maximal modulus among all the eigenvalues; (ii) a second theorem, stating 
that if a Z-matrix that is a nonsingular M-matrix is also "irreducible" , then its inverse is strictly positive. The 
irreducibility of J follows from the symmetry of the matrices Kij and Rij [9] . As the positive (Perron) eigenvector of 
— is the eigenvector of J that corresponds to its smallest magnitude eigenvalue, it defines the slowest relaxation 
mode (for a given set of temperatures). Therefore, in the evolution of temperatures given by Eqs. (4), steady-state is 
eventually approached from the zone corresponding to simultaneous temperature increments (or decrements). 

The matrix J{t) in Eqs. (5) is obtained by substituting T(o)j(t) for Tj in Eqs. (6) and (7). Then, the nonhomogeneous 
linear system with variable coefficients, Eqs. (5), can be solved by variation of parameters [9], yielding the expression: 

T(i)(i) = C/(i) /V(r)-i-F(r)rfr, (8) 
Jo 

where U{t) is a matrix formed by columns that are linearly independent solutions of the corresponding homogeneous 
equation, with the condition that U{Q) = I (the identity matrix). The difficulty in applying this formula lies in 
computing U{t), that is, in computing the solutions of the homogeneous equation. Moreover, this computation 
demands the previous computation of the solution for T(o)(t). 

Since we are only interested in the stationary solutions of the heat-balance equations rather than in transient thermal 
states, it is possible to find an expression of these solutions that is more manageable than Eq. (8). The transient 
thermal state relaxes exponentially to the stationary solution, which is a limit cycle of the nonlinear equations, 
technically speaking [8, 9]. Therefore, the stationary solution is given by the solution of Eqs. (5) with the constant 
Jacobian matrix calculated at the steady-state temperatures, which we name Tj, i = 1,...,N} This solution is 
simply [9]: 

T(i) (i) = / exp [r J] • F{t - r) dr, (9) 

with J calculated at the point T. Furthermore, the periodic stationary solution is obtained by extending the upper 
integration limit from t to infinity: 

/•oo 

T^. it) = / exp [tJ] -Fit-r) dr. (10) 
Jo 

This function is indeed periodic, unlike the one defined by Eq. (9), so it is determined by its values for t G [0,7^. 

Note that {T^.^{t)) = 0. For numerical computations, it can be convenient to express the integral from to oo as an 
integral from to T, taking advantage of the periodicity as follows: 



/ exp [r J] ■ F{t-T)dT = Y / exp [r J] • F{t -T)dT = 

oo rT 

V exp(nr J) / exp [r J] ■ F{t - t) dr = [I - exp(r J)] "M exp [r J] ■F{t-T)dT 
n=o -^0 -^0 

(the series converges because the eigenvalues of J have negative real parts). In the last integral, the argument of F 
can be transferred to the interval [0, T]: 

/ exp [rJ] • F{t -T)dT = [ exp [rJ] • F{t - r) + f exp [rJ] ■F{t-T + T) dr, 

Jo Jo Jl 

where t G Note that the one-period shift in the argument of the last F is necessary for the argument to be in 

Some remarks arc in order. First of all, we have assumed that there is one asymptotic periodic solution of the 
nonlinear Eqs. (2) and only one (a unique limit cycle). Equivalently, we have assiimcd that the perturbation series 
converges. This assumption holds in an interval of the amplitude of heat input-variations F [9]. Besides, for the 
integrals in Eq. (10) and the following equations to make sense, it is required that exp [tJ] — ^ as r ^ oo. This is 



^ The solution can also be derived as the limit of Eq. (8) in which U{t) = exp(Jt). 
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guaranteed, because the eigenvalues of J have negative real parts, as is necessary for the steady-state to be stable. 
In fact, the eigenvalues are expected to be negative real numbers and J is expected to be diagonalizable but both 
properties are not rigorously proven [9] (however, see Sect. II B). 

If J is diagonalizable, that is to say, there is a real matrix P such that P~^JP is diagonal, then the calculation of 
the integrals is best carried out on the eigenvector basis, given by the matrix P. Using this basis, Eq. (10) is expressed 
as 



rr\00 

-'(1) 



{t) =Y,Pia exp [rXa] ^ P-'Fj {t-T)dT, i = l,...,N, (11) 

« — 1 -_i 



o=l j=l 



where the first sum runs over the eigenvectors and their corresponding eigenvalues Xa- Expression (11) allows us to 
compare the contribution of the different thermal modes. In particular, for the fast modes, such that |Aa| is large, we 
can use Watson's lemma [16] to derive the asymptotic expansion: 



I 



r exp [rA„] F„{t-r)dT=^^^ + 



where Fa = Pa/Pj- When |Aa| is large, the first term suffices (unless Fa{t) is also large, for some reason); and the 
first term is small, unless Fa{t) is large. In essence, if the fast modes are not driven strongly, they can be neglected 
in the sum over a in Eq. (11). 



1. Second order perturbative equation 

For second order in e, a straightforward calculation [9] yields the following linear equation: 

T(2) = J{t) ■ r(2) + Git) , (12) 
where J{t) is the same Jacobian matrix that appears in the first-order Eq. (5) and 

^ m -^(O)J -^Wj ~ 'c' I + -^M ^(o)i ^a)i ' i = l,---,N. (13) 

The initial condition for Eq. (12) is T(-2)(0) — 0, as for Eqs. (5). Therefore, the first-order and second-order equations 
have identical solutions in terms of their respective driving terms, although G, Eqs. (13), is a known function of t 
only when the lower order equations have been solved. The integral expression, Eqs. (10), of the stationary solution 
T^-^{t) is also valid for T^-^{t), after replacing F with G and using in Eq. (13) the stationary values T'(o)(t) = T and 
T(i)(t) = T^^{t) (which make G periodic). 

It is possible to carry on the perturbation method to higher orders, and it always amounts to solving the same 
linear equation with increasingly complicated driving terms that involve the solutions of the lower order equations. 
The example of Sect. IV shows that, in a typical case, T^^{t) is a small correction to T^-^{t), and further corrections 
are not necessary. This confirms that the perturbation method is reliable for a realistic case. 

B. Heuristic linearization 

A linearization procedure frequently used in problems of radiation heat transfer [3, 6, 7] consists of using the 
algebraic identity 

T/ - T/ = (T, + Tj){Ti + T]){Ti - Tj) 
to define an effective conductance for the radiation coupling between nodes i and j. The equation 

Rij{T,^-T^)=K^:{T,-Tj), 

defines the effective conductance 
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for specified values of the node temperatures Ti and Tj. For an orbiting spacecraft, the natiiral base vahies of the 
node temperatures are the ones that correspond to the steady-state solution of the averaged equations, namely, 
Tj, i = 1, . . . ,N. \n the special case of radiation to the environment, RiTf can be replaced with linear terms Kf-Ti 
such that = ^f^Rt, for i = 1, . . . , iV. 
The resulting linear equations are: 

N 

Ci fi = Q,{t) - J2 {K^, + K^j) (T, - T,) -K^Ti, i = l,...,N, (14) 

These equations have only conduction couplings, so they are a discretization of the partial differential equations of 
heat conduction. As a linear system of ODEs, the standard form is 

T, = ^J,,r, + ^, z = l,...,iV, (15) 

j 

where J (the Jacobian matrix) is now given by: 

Jij = Cr^ {K,, + K^) , ii^^J, (16) 



N 



k=l 



(17) 



The linear system of Eqs. (15) can be solved in the standard way, yielding: 



T{t) = exp [tJ] (t{0) + ^ exp [-rJ] • ^(t) dr^ , (18) 

where we have introduced the vector q{t), with components qi{t) = Qi{t)/Ci. We can also express the solution in 
terms of the driving function F = q — (q): 

T{t) = exp [tJ] T(0) + / exp [{t - r) J] • (F(r) + {q)) dr = (19) 
Jo 

exp[tJ]T{0)+ I exp[{t-T)J]-F{T)dT + J-\exp[tJ]-I){q) (20) 
Jo 

For large t, this solution tends to the periodic stationary solution 

T°°{t)= exp[TJ]- F{t-T)dT- J-\q), (21) 
Jo 

assuming that exp [tJ] — )■ as < — )■ oo. This is a consequence of the structure of J, as in the preceding section. In the 
present case, the eigenvalues of J, beyond having negative real parts, are actually negative real numbers, as we show 
below. 

The total conductance matrix K + is symmetric but this does not imply that J is symmetric. Nevertheless, if 
we define C = diag(Ci, . . . , Cjv), the matrix C^/^ • J • C~^/^ is symmetric, because its off-diagonal matrix elements 
are: 

X -I- k'^ 

i = l,...,N, j = l,...,N, andi^j. 



Hence, the matrix C^I'^JC similar to J, has real eigenvalues. Furthermore, C^I'^JC is diagonalized by an 
orthogonal transformation; that is to say, there is an orthogonal matrix O such that 



is diagonal. Therefore, the thermal modes are actually normal; that is to say, the modes, which are the eigenvectors 
of J and hence the columns of the matrix P = C~^^^0, are related to the eigenvectors of C^/'^ JC~^^'^ , which are 
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normal and arc given by the columns of O. Alternatively, one can say that the eigenvectors of J are normal in the 
"metric" defined by C; namely, 



which can be written in matrix form as P*CP = I. Naturally, the orthogonality of modes greatly simplifies some 
computations. 

Furthermore, the symmetry of the conductance matrix implies that the sum in Eq. (14) can be written as the action 
of a graph Laplacian [17] on the temperature vector. Naturally, the graph is formed by the nodes and the linking 
conductances. A graph Laplacian is a discretization of the ordinary Laplacian and is conventionally defined with the 
sign that makes it positive semidefinite. The zero eigenvalue corresponds to a constant function, that is, a constant 
temperature, in the present case. A vector with equal components, say. equal to 1 /\fN . is the positive (Perron) 
eigenvector of the matrix. With more generality, the Laplacian of a graph can be defined as a symmetric matrix with 
off-diagonal entries that are negative if the nodes are connected and null if they are not [18] . This definition does not 
constrain the diagonal entries and, therefore, does not imply that a graph Laplacian is positive semidefinite. It can 
be made positive definite (or just semidefinite) by adding to it a multiple of the identity matrix, which does not alter 
the eigenvectors. Of course, the eigenvector corresponding to the smallest eigenvalue does not have to be constant, 
but the Perron-Frobenius theorem [15] tells us that it is positive. By this general definition of a graph Laplacian, the 
matrix — C^/^ JC~^/^ is a diff'erent Laplacian for the same graph, and Eqs. (15) contain the action of this Laplacian 
on the vector C^/'^T. Notice that this general definition of a graph Laplacian is connected with the definition of a 
Z-matrix [15] and, actually, a symmetric Z-matrix is a graph Laplacian. If such a matrix is positive definite, then it 
is equivalent to a Stieltjes matrix, namely, a symmetric nonsingular M-matrix [15]. The general Jacobian obtained 
in Sect. II A is also such that —J and also — C^/^JC~^/^ are both nonsingular M-matrices, but they need not be 
symmetric. 

To investigate the accuracy of the approximation of the radiation terms by conduction terms, let us compare the 
periodic solution given by Eq. (21) with the first-order perturbative solution found in Sect. II A, namely, T°°{t) = 
T + T^-^{t). Of course, the Jacobian matrices in the respective integrals differ, as do the temperature vectors added 

to the integrals, namely, T or —J~^{q). While T corresponds to the authentic steady-state of the nonlinear averaged 
equations, —J~^{q) corresponds to the steady-state of Eqs. (14) after averaging, which is a state without significance, 
since we have already used the set of temperatures T of the authentic steady-state to define the radiation conductances 
Kfj in Eq. (14). Therefore, the only sensible linear solution is the perturbative solution T°°{t) = T + T^^{t), even if 
we replace the Jacobian matrix given by Eq. (6) and (7) with the one given by Eq. (16) and (17). 

In our context, the notion of radiation conductance actually follows from the symmetry of the matrices CJ or 
C^/"^ JC~^/'^ . Therefore, the most natural definition of radiation conductance probably is K!f^ = 2Rij{Tf + Tj ), that 

is, the symmetrization of the term 4RijTj in Eq. (6). This symmetrization has been tested by Krishnaprakas [11], 
considering the steady-state problem for models with up to A'^ = 1237 nodes and working with various resolution 
algorithms. He found that the effect of symmetrization is not appreciable. To estimate the effect of the antisymmetric 
part of the matrix ARijTj, namely, 2Rij{Tf — T^), on the eigenvalue problem for the Jacobian, we proceed as follows. 
We formulate this eigenvalue problem in terms of the matrix C^/^ JC~^/^, so that it is an eigenvalue problem for a 
symmetric matrix perturbed by a small antisymmetric part. This problem is well conditioned, because the eigenvectors 
of the symmetric matrix (the columns of the matrix O) are orthogonal. In particular, the perturbed eigenvalues are 
still real. Furthermore, the first-order perturbation formula for the eigenvalue \a associated with an eigenvector Ca 
[16] yields: 



N 




i=l 



N 




vanishing because the perturbation matrix SA is antisymmetric. So the nonvanishing perturbative corrections begin 
at the second order in the perturbation matrix, and, in this sense, they are especially small. 
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III. FOURIER ANALYSIS OF THE PERIODIC SOLUTION 



Given that T^-^{t) is a periodic function, it can be expanded in a Fourier series. To derive this series, let us first 
introduce the Fourier series of F{t), 

oo 

F{t)= -F'M e'"™*/'^. 

m= — oo 

Inserting this series in the integral of Eq. (10) and integrating term by term, we obtain the Fourier series for T^-^(t). 
Alternatively, we can substitute the Fourier series for both T^'^^(f) and F{t) into Eqs. (5), where J is taken to be 
constant; then we can solve for the Fourier coefficients of T^-^{t). The result is 

oo 

T^^{t)= ^ e2"'"*/'^(27rim//r- J)"'-f(m), (22) 

m=— CO 



The Fourier coefficients F{m) are obtained by integration: 

F{m) = - i^(i)e-2'^'™*/'^dt. (23) 
' Jo 

Given that F{t) is a real function, 

F{-m) = F*{m). (24) 

Furthermore, {F{t)) = implies 

F(0) = 0. (25) 

So F{t) is defined by the sequence of Fourier coefficients for positive m. This sequence must fulfill the requirement 
that limm_).oo F{m) = 0, so a limited number of the initial coefficients may suffice. 

Actually, for numerical work, Eq. (23) can be conveniently replaced by the discrete Fourier transform 

n — 1 

F(m) = - ^ F{kr/n) e-^""**^/" , (26) 
" fc=o 

which only requires sampling of the values for F(t), but also only defines a finite number of independent Fourier 
coefficients, because F{m + n) = F{m). Notice that we usually have available just a sampling of the heat inputs 
at regular time intervals, rather than the analytical form of Qi{t). To calculate the exact number of independent 
Fourier coefficients provided by Eq. (26), we must take into account Eqs. (24) and (25). If n is an odd number, 
the independent Fourier coefficients F{m) are the ones with m = 1, . . . , (n — l)/2; that is to say, there are n — 1 
independent real numbers. If n is even, the independent Fourier coefficients are the ones with m = 1, . . . , n/2, and 



n — 1 

F{n/2) = -Y,{-fF{kT/n) 



is real, so there are n — 1 independent real numbers as well. For definiteness, let n be odd. Then, we can express F{t) 
as 



F{t) = 2Re 



(n-l)/2 



Of course, the values of F{t) &t t = kT/n, A; = 0, . . . , n — I, are the sampled values employed in Eq. (26), but the 
expression is valid for any t e [0, 7^ and constitutes an interpolation of the sampled values. Naturally, the higher the 
sampling frequency n, the more independent Fourier coefficients we have and the more accurate the representation of 
F{t) is. 
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As is well known, the Fourier series of a function F(t) that is piccewise smooth converges to the function, except 
at its points of discontinuity, where it converges to the arithmetic mean of the two one-sided limits [19]. However, 
the convergence is not uniform, so that partial sums oscillate about the true value of the function near each point 
of discontinuity and "overshoot" the two one-sided limits in opposite directions. This overshooting is known as the 
Gibbs phenomenon, and, in our case, produces typical errors near the discontinuities of the driving function F. These 
discontinuities are due to the sudden obstructions of the radiation on parts of the aircraft that occur at certain orbital 
positions, for example, when the Sun is eclipsed.^ Section IV B shows that the Gibbs phenomenon at eclipse points 
can be responsible for the largest part of the error of the linear method when the discrete Fourier transform is used. 

The approximation of T^-^{t) provided by the n samples of F{t) is, of course, 



r(T)(t) = 2Rc 



\n-l)/2 

J2 e'^^imt/T ^2mmI/T - J)~^ ■ F{m) 



(27) 



and is valid for any t € [0,7~]. However, if we are only interested in T^-^{t) &t t = kT/n, = 0, . . . , n — 1, we can 
compute these values with the inverse discrete Fourier transform 



r(~) [kT/n) = ^ e^"''"'^/" Urn mod (m + , v\ - I /T - j) ■ F{m) , 

m=0 V L V / J / 



(28) 



where, for m = (n -|- l)/2, . . . , n — 1, F(rn.) = F{m — n) = F*{n — m), and where mod (•,«) gives the remainder 
of the integer division by n. This inverse discrete Fourier transform can be more convenient for a fast numerical 
computation. Regarding computational convenience, the discrete Fourier transform, be it direct or inverse, is best 
performed with a fast Fourier transform (FFT) algorithm. The classic FFT algorithm requires n to be a power of 
two [20, 21]; in particular, it has to be even. 

The function T^-^{t), computed by Fourier analysis from n samples of Qi{t), is to be compared with the one 
computed by a numerical approximation of the integral formula, Eq. (10), in terms of the same samples. Naturally, 
we can use instead of the integral over r e [0,oo] the integral over re [0,T] below Eq. (10). This integral can be 
computed from the n samples of F{t) by an interpolation formula, say the trapezoidal rule. It is not easy to decide 
whether this procedure is more efficient than Fourier transforms. Considering that the substitution of the continuous 
Fourier transform, Eq. (23), by the discrete transform, Eq. (26), is equivalent to computing the former with the 
trapezoidal rule, the integral formula may seem more direct. In particular, this formula allows us to select the values 
of t for which we compute T'^'^^(t) independent of the sampling frequency, so we can choose just a few distinguished 
orbital positions and avoid the computation of all the n~ 1 integrals (one is removed by the condition {T^-^{t)) — 0). 

Note that the computation of all of the independent F{m) with Eq. (26) is equivalent to the computation of precisely 
n — 1 integrals. However, the efficiency of the FFT reduces the natural operation count of this computation, of order 

n?, to order nlogn; so its use can be advantageous, nevertheless. 

It goes without saying that the second-order perturbative contribution (t) to the stationary solution is given by 

the right-hand side of Eq. (27) with the Fourier coefficients F{m) replaced by the Fourier coefficients of the function 
G{t) defined in Sect. EAl. 



IV. TEN-NODE MODEL OF A MOON-ORBITING SATELLITE 



To test the previously explained methods, we construct a small thermal model of a simple spacecraft, namely, a 
ten-node model of a Moon-orbiting satellite. Our satellite ten-node model supports a basic thermal structure and is 
simple enough for allowing one to explicitly display the main mathematical entities, e.g., the matrices K, R and J. 
The satellite consists of a rectangular parallelepiped (a cuboid) of square base plus a small cylinder on one of its sides 
that simulates an observation instrument, as represented in Fig. 1. In addition, at a height of two thirds of the total 
height, there is an inner tray with the electronic equipment. The dimensions of the cuboid are 0.2 m x 0.2 m x 0.3 m, 
and the cylinder has a length of 0.1 m and a radius of 0.04 m. The satellite's frame is made of aluminum alloy, using 



Strictly speaking, the function Qi{t) is always continuous but it undergoes sharp variations at some times. These sharp variations can 
be considered as discontinuities, especially, if the function is sampled. 
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Node id. 6 




t: 2 mm 



FIG. 1: Satellite's structure and node description. The front face, removed to see the interior, corresponds to node 1 and is 
equivalent to nodes 3 or 4. 



Node 


G (J/K) 


m (w) 


(°C) 


1 


331.7 


15.18 


2.6 


2 


147.4 


2.30 


3.6 


3 


331.7 


15.17 


2.6 


4 


331.7 


14.80 


2.3 


5 


196.6 


3.91 


0.2 


6 


98.3 


0.63 


2.2 


7 


196.6 





6.3 


8 


31.9 


1.70 


4.7 


9 


800.0 


4.35 


15.9 


10 


1400.0 


6.15 


11.1 



TABLE I: Node capacities and mean heat inputs with their associated steady-state temperatures. 



plates 1 mm thick, except the bottom plate, which is 2 mm thick. This plate plays the role of a radiator and its outer 
surface is painted white to have high solar reflectance. The cylinder is made of the same aluminum alloy, as well as 
the tray; they are 0.5 mm and 2 mm thick, respectively. The sides of the satellite, except the one with the instrument, 
are covered with solar cells, which increase the sides' thickness to 2.25 mm. 

The thermal model of the satellite assigns one node to each face of the cuboid, one more to the cylinder and 
another to the tray, that is, eight nodes altogether. Furthermore, to conveniently split the total heat capacitance of 
the electronic equipment, it is convenient to add two extra nodes with (large) heat capacitance but with no surface that 
could exchange heat by radiation. Nodes of this type are called "non-geometrical nodes" . In the present case, they 
represent two boxes with equipment placed above and below the tray, respectively. We order the ten nodes as shown 
in Fig. 1. The lower box (node 10) is connected to the radiator by a thermal strap. Given the satellite's structure and 
assuming appropriate values of the specific heat capacities, it is possible to compute the capacitances C'i, i — 1 . . . , 10, 
with the result given in Table I. Using the value of the aluminum- alloy heat conductivity and assuming perfect contact 
between plates, we compute the conduction coupling constants Kij between nodes i,j = 1,...,8. The remaining 
conduction coupling constants are given reasonable values, shown in Eq. (29). The computation of the radiation 
coupling constants Rij, i,j = 1, . . . , 8, and Ri, i = 1, . . . , 8, and indeed the computation of the external radiation 
heat inputs requires a detailed radiative model of the satellite, consisting of the geometrical view factors and the 
detailed thermo-optical properties of all surfaces. This radiative model allows us to compute the respective absorption 
factors [1]. 

The thermo-optical properties of the surfaces are assumed to be as realistic as possible, given the simplicity of the 
thermal model. All radiation reflection is assumed to be diffuse, as is common for many types of surfaces. The inner 
surfaces are painted black and have high emissivity, e = 0.84, to favor the uniformization of the interior temperature. 
The outer surfaces are of three types. The three sides covered with solar cells also have high emissivity, e = 0.84, to 
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FIG. 2: Thermo-optical properties of the satellite's surfaces (details are given in the text). 



favor the cooling of the solar cells. On the other hand, they have high solar absorptivity, as = 0.75. Of this 0.75, 0.18 is 
processed into electricity and the remaining 0.57 dissipates as heat in the solar cells. The top surface, the surface with 
the cylinder, and the cylinder itself (its two sides) have low emissivity, e = 0.1, and low solar absorptivity, ctg = 0.2, 
which are chosen to simulate the effect of a multilayer insulator. In contrast, the bottom surface simulates a radiator, 
with e — 0.8 and as — 0.2 (like an optical solar reflector). All of these thermo-optical properties are summarized 
in Fig. 2. For the computation of the corresponding absorption factors, we employ the ray-tracing Monte-Carlo 
simulation method provided by ESARAD"^^ (ES ATAN'^'^ 's radiation module) [14]. 

Taking into account the above information, one obtains the following conduction (in W/K) and radiation (in W/K*) 
matrices: 



1 

10 



/ 3.47 5.64 2.86 2.00 4.50 

3.47 3.47 1.67 1.33 3.50 3.00 

3.47 5.64 2.86 2.00 4.50 

5.64 5.64 2.86 2.00 4.50 

2.86 1.67 2.86 2.86 3.00 

2.00 1.33 2.00 2.00 

4.50 3.50 4.50 4.50 4.50 6.00 

3.00 00000000 

4.50 

V 3.00 6.00 

/ 5.06 4.63 5.05 3.68 2.71 6.39 \ 
5.06 5.05 4.63 3.68 2.70 6.39 0.13 



(i?.,) = 10 



-10 



4.63 5.05 5.06 3.69 2.71 6.39 
5.05 4.63 5.06 3.69 2.70 6.38 
3.68 3.68 3.69 3.69 3.57 
2.71 2.70 2.71 2.70 7.19 
6.39 6.39 6.39 6.38 3.57 7.19 
0.13 
00000000 
00000000 

















0/ 



(Ri) = lO"'' (2.86, 0.32, 2.86, 2.86, 1.81, 0.23, 0,0.23, 0,0). 



(29) 



(30) 



(31) 



The satellite's thermal characteristics are defined by the data set {Ci, Kij, Rij, Ri}, but the radiation heat exchange 
depends on the nodal temperatures, which in turn depend on the heat input. As explained in Sect. II A, the appropriate 
set of nodal temperatures corresponds to the steady-state for averaged heat inputs, given by the algebraic equation 
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FIG. 3: The 111 positions of tlie satellite in its orbit. The sunlight comes along the x axis. 

that results from making r(o)i = in Eq. (4). Since we need the external heat inputs and, therefore, the orbit, we 
proceed to define the orbit characteristics. 

We choose a circular equatorial orbit 26926 m above the Moon's surface, such that T = 6660 s. The radiation 
heat input to the satellite consists, on the one hand, of the solar irradiation and the Moon's albedo, and, on the 
other hand, of the Moon's constant IR radiation. We take 0.12 for the mean Moon's albedo and 270 K for the black- 
body equivalent temperature of the Moon. There is also heat produced by the dissipation of electrical power in the 
equipment (nodes 9 and 10). For the sake of simplicity, the dissipation rate is assumed to be constant, equal to the 
mean electrical power generated in an orbit. In a part of the orbit, the Moon eclipses the Sun, so the satellite receives 
no direct sunlight or albedo, although there is always IR radiation from the Moon. The satellite is stabilized such 
that the cylinder (the "observation instrument" ) always points to the Moon and the longer edges are perpendicular 
to the orbit. The radiation heat input can be computed by taking into account the given orbital characteristics 
and the satellite's thermo-optical characteristics, in particular, the absorption factors. It has been computed with 
ESARAD"""^, taking 111 positions on the orbit, that is, at intervals of one minute. 

In Fig. 3, all 111 positions are plotted. The initial position of the satellite is at the subsolar point and it moves 
towards the east. The total external radiation heat input to the first eight nodes (the ones that receive radiation) is 
plotted in Fig. 4 (only at every other position, for clarity). Note the symmetry between nodes 1 and 3, which denote 
the lateral faces, covered with solar cells. Node 4 corresponds to the back side, also covered with solar cells. So the 
external radiation load on it has a similar time variation, but it is displaced. The solar radiation absorbed by all solar 
cells results in an orbital mean power rate of 10.5 W, dissipated in the equipment and split between nodes 9 and 10, 
which receive 4.35 and 6.15 W, respectively. The external radiation absorbed by the side with the cylinder (node 2) is 
considerably smaller than the radiation absorbed by the sides with solar cells, due to the low value of as (and of e, as 
well) for the corresponding surface. The bottom and top outer surfaces, which belong to nodes 5 and 6, respectively, 
have view factors for the external radiation that are much less favorable than those of the side surfaces. Nevertheless, 
the amount of lunar IR radiation absorbed by the bottom surface, due to its high e, is such that the orbital mean of 
the external heat input to node 5 is, in fact, larger than the one for node 2 (see Table I). Naturally, node 7, with no 
outer surfaces, does not absorb any external radiation. 

To determine the hot and cold cases of the orbit, we compute the total heat load on the satellite for each position 
in the orbit, finding a maximum of 90.59 W at position 14 and a minimum of 18.65 W at any position in the eclipse, 
during which all the heat loads stay constant. The solution of the corresponding steady-state problems at position 14 
and at a position in the eclipse yields the two sets of nodal temperatures (for the given node order): 

Hot {29.0, 35.5, 41.2, 38.5, 30.2, 35.1, 39.1, 35.0, 48.7, 42.9} °C. 

Cold {-46.8, -45.0, -46.8, -49.1, -45.7, -47.1, -42.8, -44.1, -33.2, -37.0} °C. 

The results in Sect. IV A show that the periodic thermal state does not reach these extreme temperatures, which 
could endanger the performance of the satellite. 




FIG. 4: Variation of the external heat input with time t (in minutes) along the orbit. For both plots, the node numbers are 
denoted by shape, in the order: dots (1,5), squares (2,6), diamonds (3,7), and triangles (4,8). 



A. Jacobian matrix and periodic solution 



We compute the averages (Qi) and substitute for them and the data set {Kij,Rij, Ri} in Eq. (4) to find the steady 
state temperatures. The results are given in Table I. Then, according to Eqs. (6) and (7), the Jacobian matrix is (in 



J = 10-''' 



-6.99 


1.18 


0.12 


1.83 


0.95 


0.67 


1.52 











2.64 


-12.93 


2.64 


0.26 


1.33 


1.06 


2.75 


2.04 








0.12 


1.17 


-6.99 


1.83 
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0.67 


1.52 











1.83 


0.12 


1.83 


-7.64 


0.95 


0.67 


1.52 











1.61 


1.01 


1.61 


1.61 


-8.26 





0.16 








1.53 


2.27 


1.59 


2.27 


2.27 





-9.20 


0.64 











2.56 


2.06 


2.56 


2.56 


0.15 


0.31 


-15.60 





2.29 


3.05 





9.43 

















-10.05 























0.56 





-0.56 

















0.21 





0.43 








-0.64 



By inspection, one can check that it has nonnegative off-diagonal and negative diagonal elements, that is to say, — J 
is a Z- matrix. It is also diagonally dominant, namely, \Jii\ > '^j^i \Jij\- The eigenvalues of J are 

-10"^{182.20, 154.30, 103.40, 98.03, 86.12, 71.09, 71.04, 14.90, 5.70, 1.72}. 

Their inverses (in absolute value) give us the typical relaxation times of the corresponding thermal modes. Thus, we 
deduce that relaxation time of the fastest mode is about 55 s, whereas the relaxation time of the slowest one is 5813 
s. The latter time is similar to T = 6660 s. 

The eigenvalues are real numbers and, furthermore, J is diagonalizable, because the eigenvalues are different. 
Both properties also follow from C^/^JC~^^^ being almost symmetric: its antisymmetric part, SA = {C^/^JC~^^^ — 
C-i/V*Ci/2)/2, is relatively smah, namely, \\6A\\/\\C^/'^JC"'^/'^\\ < 10~^, where the matrix norm is the Frobenius 
norm (other standard matrix norms yield similar values). Therefore, the notion of "radiation conductance" (Sect. II B) 
is appropriate in this case, as concerns its use in the linear equations. The thermal modes are almost normal, namely, 
the eigenvector matrix P is such that P^CP — I with an error < 0.002. The most interesting eigenvector of J is, of 
course, the positive (Perron) eigenvector, which corresponds to the slowest mode. The normalized positive eigenvector 
is 

(0.259, 0.276, 0.259, 0.257, 0.275, 0.267, 0.327, 0.264, 0.471, 0.423) K. 



Note that the temperature increments are of a similar magnitude, except the ones of node 7 and, especially, nodes 
9 and 10, which are associated, respectively, for the tray and the boxes of electronic equipment. The next mode, 
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7;°° (°C) T,^ (°C) 




t 



(a) nodes 1-5 (b) nodes 6-10 



FIG. 5: Variation of the ten nodal temperatures with time t (in minutes). For both plots, the node order is: dots (1,6), squares 
(2,7), diamonds (3,8), upward triangles (4,9), and downward triangles (5,10). 

corresponding to the eigenvalue —5.70 • 10~^, has one negative component (the ninth), and the remaining modes have 
more than one. 

To calculate T^-^{t), we choose the Fourier series of Eq. (27) or, rather, the inverse discrete Fourier transform of 

Eq. (28), which can be computed with a FFT algorithm. The Fourier coefficients F{m) can also be computed with 
the FFT, according to Eq. (26). Once the vector at the 111 positions is available, the set of nodal temperatures 

corresponding to the first-order perturbative solution is T°°{t) = Ti + T^^-{t), i = 1, . . . , 10, plotted in Fig. 5. A 
measure of the accuracy of this perturbative calculation is given by the second-order calculation in the next section. 
The truncation of the Fourier series imposed by the sampling of F also is a source of error, unrelated to perturbation 
theory. The piecewise smoothness of the function T^-^{t) suggests that the error is small (but see Sect. IV B). 

It is also interesting to see if the first-order perturbative calculation is affected by neglecting the fastest modes: 
according to the analysis at the end of Sect. II A, these modes are expected to contribute in proportion to their 
relaxation times. The fastest mode relaxes in about 55 s, a short but non-negligible time. As a consequence, its 
contribution to T^^y which we find to have a maximum magnitude of 0.8 K, is small but non-negligible. But we can 
deduce that still faster modes, which would appear in a thermal model of the satellite with more nodes, are hardly 
necessary. 

From the engineering standpoint, note that this satellite thermal model is successful, insofar as it predicts that 
all nodal temperatures stay within adequate ranges. In particular, nodes 9 and 10, corresponding to the boxes with 
electronic equipment, stay within the range from 4 to 23 °C. These nodes are inner nodes with large thermal capacity 
and, hence, are protected against the larger changes in the external heat inputs. In contrast, the outer nodes are very 
exposed and undergo considerable variation in temperature, with especially sharp changes at the beginning and end 
of the eclipse. 

1. Second-order correction 

According to Sect. II A 1, the second-order perturbative correction T^-^ to the periodic stationary solution is obtained 

by the same procedure as that for , but using a different driving function G that is computed from T and from 
T^-^ itself. The computations are straightforward and they yield the correction plotted in Fig. 6. This correction is 
always negative, because the negative term in the expression for G, Eq. (13), dominates over the positive term. The 
equation (10) for T^-^ and the corresponding equation for T^-^ are both linear, so T^-^ and are proportional to 
the respective driving functions; and we can compare their magnitudes by comparing those driving functions, say, 
comparing typical values of Q and QRiff T^-^^ . This latter quantity can be roughly estimated as 6 • 2 • 10~^ • 300^ • 20^ 
W ~ 0.4 W, whereas Q ~ 10 W (Table I). Their ratio is about 25, which roughly agrees with the ratio of T^^^ to T^y 
as can be seen by comparing Fig. 5 to Fig. 6. 

The order of magnitude of the second-order correction suggests that higher orders are not necessary, as we show 
next. 
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FIG. 6: Second-order correction to the temperatures (same node order as in Fig. 5). 



B. Direct integration of the nonlinear equations 



Of course, the periodic solution T°°{t) can also be obtained by direct integration of the nonlinear Eqs. (1) with an 
adequate solver, based on Runge-Kutta or other methods [13]. Since the nonlinear analysis proves that the solution 
of Eqs. (1) will converge to the stationary periodic solution [9], the numerical solver must be run until this periodicity 
is established. Periodicity can be enforced by comparing the nodal temperatures at the beginning and the end of 
every period and demanding that they be equal, within some tolerance. To do this, ESATAN"'"'^ provides the routine 
SOLCYC [14]. We have employed this routine to obtain the nonlinear equations' cyclic solution (which is established 
in 10 periods if the tolerance is set to 0.001). 




FIG. 7: The error in the temperature computed up to the second order. Only showing nodes 1 (dots), 3 (squares), and 8 
(diamonds), which give rise to the largest errors. 

To compare the cyclic solution obtained by the linear method with the "exact" solution T°°{t) obtained by SOLCYC, 
we quantify the deviation by the vector of "errors" 

AT{t) = T°°(<) -\f + T[^^{t) + T^~^{t) 

The largest components of AT are plotted in Fig. 7 (the remaining components stay in the range [—0.1, 0.1] °C). 
There seem to be two branches for each node, but, in fact, it is an effect produced by the high-frequency oscillations 
of each ATi{t). Notice that the error is generally small when compared with T^^{t); namely, [Ar(t)[ ^ [T(5^^(t)[ for 
each t = kT/n, fc — 0, . . . , n — 1, except near the two positions corresponding to the beginning and end of the eclipse, 
where the heat inputs have discontinuities (Fig. 4). These discontinuities induce oscillations of certain amplitude 
about the true values of T°°(t), due to the Gibbs phenomenon in the discrete Fourier transform (Sect. III). To 
suppress the oscillations, we would need a specific method for treating the Gibbs phenomenon; for example, we could 
use a smoother method for Fourier series summation, such as Cesaro summation [19] 
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V. SCALABILITY AND COMPLEXITY OF THE LINEAR METHOD 

The ten-node model studied in Sect. IV is too small to pose a computational problem, whether we employ the 

linear method or directly integrate nonlinear Eqs. (1). To assess the practical applicability of the linear method, we 
must study how it scales to realistic sizes and then compare it with the direct integration of Eqs. (1). Naturally, 
the first step in the method is to compute the steady-state temperature Tj , but this computation, arguably, is not 
a substantial part of the whole process and we do not consider it (it may take, say, between a few percent and one 
fith of the whole process, depending on the circumstances). The size of an orbiting spacecraft thermal model can be 
scaled with respect to the number of nodes, N, or the number of different positions in the orbit, n. However, while 
realistic models must have many nodes, an n of about one hundred can generally be suitable. Note, in particular, 
that the condition that ^ n in the ten-node model is likely to be reversed in realistic models. For these reasons, 
the scaling with respect to N is more relevant. 

For the moment, let us neglect any special fcatin-c of Eqs. (1), such as possible coefficient sparsity. Then, the 
computational complexity of numerically integrating those equations is of order N'^k, where k is the nimiber of time 
steps taken. On the other hand, the complexity of the numerical integration of linear Eqs. (5) is also of order N'^k 
[excluding the computation of J{t)]. We can employ the explicit integral, Eq. (8), which can be calculated with 
the trapezoidal rule, for example, but this does not reduce the complexity of the computation. However, for the 
stationary solution, Eq. (10), its expression as an integral over a period T sets the number of time steps to n, which 
is advantageous if n < k, where k now is the number of steps necessary for some initial temperatures to relax to the 
stationary solution. Since we have to evaluate the integral in Eq. (10) for several times t, it is preferable to use the 
FFT, as discussed in Sect. Ill, so that the total operation count is of order N"^ nlogn. 

However, we have not taken into account matrix operations of considerable complexity. For example, the discrete 
Fourier transform, Eq. (27), involves the inverse of an N x N matrix (—J plus a multiple of the identity matrix), 
and the inversion of a matrix is generally a process of order A'^'^ [20]. If we have to employ a process of order N^, 
we may as well diagonalize J, because the diagonalization of a matrix also is generally a process of order [20] 
and the diagonalization of J has several uses. Let us assume that we carry out this diagonalization and determine 
the independent thermal modes, which can then be employed to express Eq. (10) as Eq. (11) or to simplify the 
matrix operations in Eq. (27). If we roughly compare the computational complexity of order A''^ with the numerical 
integration complexity of order N'^k, we deduce that the diagonalization is worthwhile when N < k. Using the above 
estimate n ~ 100 and taking as relaxation time k ^ 5n (the rough value for the ten- node model), we deduce that 
the determination of thermal modes can be useful just as a computational procedure for models with a few hundred 
nodes. 

Let us now consider that J is surely a sparse matrix, so iterative matrix methods can take advantage of this 
characteristic. In fact, there are two degrees of sparsity in J, associated with conduction or radiation coupling terms. 
The conductance matrix K has to be very sparse, because conduction is a local process, so each node can only 
be coupled to a few nodes. In contrast, radiation is a nonlocal process and couples any pair of nodes that has a 
nonvanishing view factor. In addition to the different sparsity of conduction and radiation coupling matrices, there 
are two other circumstances that make them different: (i) we are assuming that the conductive coupling matrix just 
depends on material properties, so it is independent of the reference temperatures Tf, (ii) the conduction coupling 
terms are significantly larger than the radiation coupling terms for the natural values of those temperatures. All of 
this suggests separating the conduction and radiation parts of J, and then diagonalizing the conduction part. 

Based on the study in Sect. II B, the diagonalization of the conduction part of J boils down to the diagonalization 
of the (generalized) Laplacian matrix —C~^^'^KC~^^'^ , and this matrix is sparse. Suitable iterative algorithms to 
perform this diagonalization are, for example, the Lanczos [21] or the Davidson [22] algorithms. These algorithms are 
particularly useful when only a few of the largest or smallest eigenvalues are needed. This is indeed our case, as only 
the slower modes are expected to contribute; to T^-^ and T^^^ . Once the conduction part of J has been diagonalized, in 
the sense that the lowest eigenvalues and eigenvectors of the corresponding Laplacian matrix are known, the radiation 
part of J can be treated as a perturbation, using matrix perturbation methods [16, 21]. 

Iterative matrix methods, combined with matrix perturbation methods or other methods, if necessary, can reduce 
the order N'^ to iV^ or even to almost linear and so allow us to diagonalize the Jacobians for the largest values of N 
that appear in current thermal spacecraft models. Of course, the sparsity of thermal coupling matrices also facilitates 
the direct integration of the nonlinear Eqs. (1). However, their nonlinearity prevents one from taking advantage 
of the above mentioned approximations methods, for example, the reduction to the small subspace of slow modes, 
or the splitting into conduction and radiation in which the latter is treated as a matrix perturbation. Moreover, 
the linearization is useful in various respects. For example, the overall relaxation time, given by the eigenvalue of 
smallest magnitude, can be effectively bounded by inequalities [23] and some of these bounds can be found with little 
computational effort. 
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VI. SUMMARY AND DISCUSSION 

We have studied the evolution of the thermal state of an orbiting spacecraft and developed a linear approach 
to this problem that is based on a rigorous pcrturbative treatment of the exact nonlinear equations. The first- 
order perturbation equations, Eqs. (5), constitute the basic linear system, which can be applied to higher orders 
after calculating the corresponding driving terms. As the Jacobian matrix of the nonlinear equations has negative 
eigenvalues, the linear equations describe the relaxation to a stationary thermal state, namely, a periodic solution 
that is independent of the initial conditions and only depends on the external heat input. This relaxation is similar 
to the relaxation to steady-state under constant external heat load. 

We have shown that the perturbative treatment reveals the scope of a common linearization procedure of a heuristic 
nature, in which the nonlinear equations are rendered linear by the definition of radiation conductances (Sect. II B). If 
one previously calculates, with the correct nonlinear Eqs. (4), the reference steady-state condition that corresponds to 
the average external heat input, the deviation from that steady-state is well approximated by the linear equations with 
radiation conductances. The Jacobian matrix corresponding to radiation conductances, obtained in Eqs. (16) and 
(17), is related to a symmetric matrix and, therefore, is easier to diagonalize. Furthermore, this relation implies that 
the thermal modes are normal, like the vibrational modes of a mechanical system. Although the notion of radiation 
conductance is just an approximation, it serves nonetheless to show that the Jacobian matrix is diagonalizable and 
has real eigenvalues. 

The diagonalization of the Jacobian matrix is useful for the computation of the stationary thermal state and also 
provides information on the relaxation to that state, because the relaxation times of the thermal modes are the 
inverses of the eigenvalues. These times span a considerable range, but the longest times are much more significant 
than the shortest times, because the latter depend on the details of the lumped-parameter thermal model employed 
whereas the former are essentially independent of it. In fact, a thermal model that has more nodes and therefore 
more details also has more thermal modes; but the slowest modes, which correspond to temperature changes in large 
parts of the spacecraft, are hardly affected by the details, whereas the fast modes can be significantly altered. The 
slowest mode, in particular, corresponds to a simultaneous but non-uniform increase (or decrease) of the temperature 
throughout the spacecraft and is hardly altered by small-scale changes. 

The computation of the stationary thermal state with the linear method relies on an explicit integral, Eq. (10), or 
a Fourier expansion, Eq. (27). Given a sampling of the thermal driving function at equal time intervals, the periodic 
solution can be obtained through two discrete Fourier transforms: a direct transform to get the Fourier coefficients of 
the driving function and an inverse transform of the coefficient vector multiplied by a suitable matrix (Sect. III). Of 
course, the discrete Fourier transforms are best performed with a EFT algorithm. This computation is more efficient 
than the numerical computation of the integral, Eq. (10), if we need the values of the temperatures at all the given 
sampling times. However, the Fourier transform presents the Gibbs phenomenon, associated with sudden variations 
of the heat loads, as occur at eclipse times, for example. The Gibbs phenomenon introduces errors, but these errors 
could be suppressed with special methods. 

The computation of the thermal modes and the stationary thermal state for a satellite ten-node thermal model 
confirms the validity of the linear method for a minimal but realistic model. The relaxation times span a considerable 
range, between 55 seconds and nearly one hundred minutes. Of course, the latter time must be almost independent 
of the particular thermal model used, whereas the former has no intrinsic significance, and, if the number of nodes 
grew, that time would shrink (thus further expanding the range of relaxation times). The slowest mode corresponds 
to node temperature increments with the same sign (positive by convention), whereas the increments corresponding 
to other modes have both signs. The periodic variation in the external heat input (Fig. 4) excites the thermal modes 
and produces a definite pattern of stationary temperature oscillations, well approximated by the first-order solution 
(Fig. 5). The second-order correction is small compared to the first-order solution, but it is worth computing, as it 
reaches 1.7 K. Higher order corrections are essentially negligible, but the error due to the Gibbs phenomenon at the 
eclipse positions reaches 0.6 K (at the most). 

Focusing on the computational aspects of the linear approach, we have studied how it scales with the number N 
of nodes and the number n of sampling positions on the orbit. If the Jacobian matrix is dense, the complexity of the 
corresponding matrix operations is of order N^. It is convenient to employ just one matrix operation, namely, the 
diagonalization of the Jacobian matrix, because then only a few of the slowest modes are needed for the remaining 
operations, so these have negligible complexity. The complexity of a direct numerical integration of the nonlinear 
equations is of order N'^k, k being the number of time steps necessary for relaxation. For a low-altitude orbit, k is 
expected to be on the order of one thousand, as for our Moon-orbiting satellite. Therefore, the linear method would 
be computationally effective as just an integration method only for models with a few hundred nodes. At any rate, 
the Jacobian matrix can be assumed to be sparse, and its conduction part can be assumed to be especially sparse, in 
addition to being the larger part of the Jacobian matrix and also being independent of the orbit. As the orbit may 
be subjected to changes in the planning of a mission, a convenient strategy probably is to diagonalize the conduction 



19 



part at the outset and, when needed, add the radiation part within some approximation scheme. This strategy can 
be far more efficient than integrating the nonlinear equations each time. 

Moreover, the strcngtli of the linear approach lies with the insight that it provides about the thermal behavior of 
the spacecraft, as embodied by the decomposition of its thermal modes, of which only the slowest ones are significant. 
These significant modes can actually be obtained with a reduced thermal model using few nodes. Therefore, the linear 
approach is especially useful in the context of reduced models. Furthermore, it provides a method for model reduction 
based on the mode decomposition: this decomposition can be used to group nodes. Indeed, there is a technique for 
graph partitioning based on the eigenvalues and eigenvectors of the Laplacian matrix of the graph [18, 24]. According 
to Sect. II B, this technique is applicable to the Jacobian matrix, but the details of this application are beyond the 
scope of the present paper and are left for future work. 

Finally, our linear approach can surely be applied to other cyclic heating processes that involve radiation heat 
transfer. 
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